Mouse analysis enge
Jason T. Serviss
2019-10-02
packages <- c("CIMseq", "CIMseq.testing", "CIMseq.data", "tidyverse", "ggthemes")
purrr::walk(packages, library, character.only = TRUE)
rm(packages)
##DATA
load('../data/CIMseqData.rda')
load('../data/sObj.rda')
if(!dir.exists('../figures')) dir.create('../figures')p <- plotUnsupervisedClass(cObjSng, cObjMul, palette('total'))
p
ggsave(
plot = p,
filename = '../figures/MGA.enge.classes.pdf',
device = cairo_pdf,
height = 180,
width = 180,
units = "mm"
)
p <- plotUnsupervisedMarkers(
cObjSng, cObjMul,
c("Lgr5", "Ptprc", "Chga", "Dclk1", "Alpi", "Slc26a3", "Atoh1", "Lyz1"),
pal = RColorBrewer::brewer.pal(8, "Set1")
)
p
ggsave(
plot = p,
filename = '../figures/MGA.enge.markers.pdf',
device = cairo_pdf,
height = 180,
width = 180,
units = "mm"
)
p <- plotUnsupervisedMarkers(cObjSng, cObjMul, "Mki67")
p
ggsave(
plot = p,
filename = '../figures/MGA.enge.Mki67.pdf',
device = cairo_pdf,
height = 180,
width = 180,
units = "mm"
)
# p <- plotUnsupervisedMarkers(cObjSng, cObjMul, "Hoxb13")
# ggsave(
# plot = p,
# filename = '../figures/MGA.enge.Hoxb13.pdf',
# device = cairo_pdf,
# height = 180,
# width = 180,
# units = "mm"
# )Mouse age
p <- getData(cObjSng, "dim.red") %>%
matrix_to_tibble("sample") %>%
inner_join(MGA.Meta, by = "sample") %>%
mutate(subject_age = as.character(subject_age)) %>%
mutate(subject_age = parse_factor(subject_age, levels = sort(unique(subject_age)))) %>%
ggplot() +
geom_point(aes(V1, V2, colour = subject_age), alpha = 0.7) +
ggthemes::theme_few() +
labs(x = "UMAP dim. 1", y = "UMAP dim. 2") +
scale_colour_manual(values = rev(col40())) +
guides(colour = guide_legend(title = "Age (days)"))
p
ggsave(
plot = p,
filename = '../figures/MGA.enge.QC.age.pdf',
device = cairo_pdf,
height = 180,
width = 180,
units = "mm"
)Mouse sex
p <- getData(cObjSng, "dim.red") %>%
matrix_to_tibble("sample") %>%
inner_join(MGA.Meta, by = "sample") %>%
mutate(subject_sex = as.character(subject_sex)) %>%
mutate(subject_sex = parse_factor(subject_sex, levels = sort(unique(subject_sex)))) %>%
ggplot() +
geom_point(aes(V1, V2, colour = subject_sex), alpha = 0.7) +
ggthemes::theme_few() +
labs(x = "UMAP dim. 1", y = "UMAP dim. 2") +
scale_colour_manual(values = rev(col40())) +
guides(colour = guide_legend(title = "Sex"))
p
ggsave(
plot = p,
filename = '../figures/MGA.enge.QC.sex.pdf',
device = cairo_pdf,
height = 180,
width = 180,
units = "mm"
)Total counts
cs <- colSums(getData(cObjSng, "counts"))
p <- getData(cObjSng, "dim.red") %>%
matrix_to_tibble("sample") %>%
inner_join(tibble(sample = names(cs), counts = cs)) %>%
ggplot() +
geom_point(aes(V1, V2, colour = counts), alpha = 0.7) +
ggthemes::theme_few() +
labs(x = "UMAP dim. 1", y = "UMAP dim. 2") +
scale_colour_viridis_c() +
guides(colour = guide_colourbar(title = "Total counts"))## Joining, by = "sample"p
ggsave(
plot = p,
filename = '../figures/MGA.enge.QC.counts.pdf',
device = cairo_pdf,
height = 180,
width = 180,
units = "mm"
)Total genes
tg <- apply(getData(cObjSng, "counts"), 2, function(c) length(which(c != 0)))
p <- getData(cObjSng, "dim.red") %>%
matrix_to_tibble("sample") %>%
inner_join(tibble(sample = names(tg), genes = tg)) %>%
ggplot() +
geom_point(aes(V1, V2, colour = genes), alpha = 0.7) +
ggthemes::theme_few() +
labs(x = "UMAP dim. 1", y = "UMAP dim. 2") +
scale_colour_viridis_c() +
guides(colour = guide_colourbar(title = "Total genes"))## Joining, by = "sample"p
ggsave(
plot = p,
filename = '../figures/MGA.enge.QC.genes.pdf',
device = cairo_pdf,
height = 180,
width = 180,
units = "mm"
)Cell type markers
markers <- c(
"Lgr5", "Ptprc", "Chga", "Dclk1", "Alpi", "Slc26a3", "Atoh1", "Lyz1", "Mki67",
"Hoxb13", "Plet1", "Osr2", "Reg4", "Sval1"
)
gene.order <- c(
"Hoxb13", "Osr2", "Atoh1", "Reg4", "Sval1", "Plet1", "Mki67", "Lgr5", "Alpi",
"Slc26a3", "Lyz1", "Dclk1", "Chga", "Ptprc"
)
d <- getData(cObjSng, "counts.cpm")[markers, ] %>%
t() %>%
as.data.frame() %>%
rownames_to_column("Sample") %>%
as_tibble() %>%
gather(gene, cpm, -Sample) %>%
inner_join(tibble(
Sample = colnames(getData(cObjSng, "counts")),
Classification = getData(cObjSng, "classification")
), by = "Sample") %>%
mutate(Classification = parse_factor(Classification, levels = classOrder.MGA('total'))) %>%
mutate(gene = parse_factor(gene, levels = gene.order)) %>%
arrange(Classification, gene) %>%
mutate(Sample = parse_factor(Sample, levels = unique(Sample))) %>%
mutate(type = case_when(
str_detect(Classification, "^S") ~ "Small intestine",
str_detect(Classification, "^C") ~ "Colon",
TRUE ~ "Other"
)) %>%
mutate(type = parse_factor(
type,
levels = c("Small intestine", "Colon", "Other")
))
under <- d %>%
mutate(idx = 1:nrow(.)) %>%
group_by(Classification) %>%
summarize(
min = (min(idx)), max = (max(idx) - 1), median = median(idx)
) %>%
ggplot() +
geom_segment(
aes(x = min, xend = max, y = 0, yend = 0, colour = Classification)
) +
geom_text(
aes(median, y = 0, label = Classification),
angle = 90, nudge_y = -0.01, hjust = 1, colour = "grey10"
) +
scale_colour_manual(values = palette('total')[classOrder.MGA('total')]) +
ylim(c(-1, 0)) +
theme_few() +
scale_x_continuous(expand = c(0, 0)) +
theme(
axis.text = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank(),
panel.border = element_blank()
) +
guides(colour = FALSE)
over <- d %>%
ggplot() +
geom_bar(aes(Sample, cpm, fill = Classification), stat = "identity") +
facet_grid(gene ~ type, scales = "free", space = "free_x") +
scale_fill_manual(values = palette('total')[classOrder.MGA('total')]) +
theme_few() +
theme(
axis.text.x = element_blank(),
axis.title.x = element_blank(),
axis.ticks.x = element_blank(),
strip.background.x = element_rect(fill = "white")
) +
labs(y = "CPM") +
guides(fill = FALSE)
library(patchwork)
p <- over + under + plot_layout(ncol = 1, heights = c(2, 1))
p
ggsave(
plot = p,
filename = '../figures/MGA.enge.geneBar.pdf',
device = cairo_pdf,
height = 300,
width = 300,
units = "mm"
)Stemness dotplot
markers <- c(
"Wnt3", "Wnt2b", "Wnt4", "Wnt5a", "Lgr5", "Axin2", "Ascl2", "Egf", "Notum",
"Kit", "Dll1", "Dll4", "Reg4", "Plet1"
)
gene.order <- markers
#scale.func <- scale_size
scale.func <- scale_radius
scale.min <- NA
scale.max <- NA
cOrderAlt <- c(
"SI Goblet", "SI Paneth", "SI Stem", "SI Transit amplifying",
"SI Progenitor early", "SI Progenitor late", "SI Enterocytes",
"C Goblet proliferating", "C Goblet Plet1 1", "C Goblet Plet1 2", "C Stem 1",
"C Stem 2", "C Stem 3", "C Transit amplifying", "C Progenitor",
"C Colonocytes", "C Goblet Plet1 neg.", "Enteroendocrine", "Tufft", "Blood"
)
p <- getData(cObjSng, "counts.cpm")[markers, ] %>%
t() %>%
as.data.frame() %>%
rownames_to_column("Sample") %>%
as_tibble() %>%
gather(gene, cpm, -Sample) %>%
inner_join(tibble(
Sample = colnames(getData(cObjSng, "counts")),
Classification = getData(cObjSng, "classification")
)) %>%
mutate(type = case_when(
str_detect(Classification, "^S") ~ "Small intestine",
str_detect(Classification, "^C") ~ "Colon",
TRUE ~ "Other"
)) %>%
group_by(gene, Classification) %>%
summarize(mean = mean(cpm), pct = 100 * (length(which(cpm != 0)) / n())) %>%
mutate(scaled.mean.exp = scale(mean)) %>%
ungroup() %>%
mutate(type = case_when(
str_detect(gene, "^Wnt") ~ "Wnt ligand",
gene %in% c("Lgr5", "Axin2", "Ascl2") ~ "Wnt target",
gene %in% c("Egf", "Notum", "Dll1", "Dll4") ~ "Other stemness factors",
gene %in% c("Plet1", "Reg4", "Kit") ~ "Stem-adjacent\nsecretory cell factors",
TRUE ~ "NA"
)) %>%
mutate(tissue = case_when(
str_detect(Classification, "^SI") ~ "Small intestine",
str_detect(Classification, "^C") ~ "Colon",
TRUE ~ "Other"
)) %>%
mutate(tissue = parse_factor(tissue, levels = c("Small intestine", "Colon", "Other"))) %>%
mutate(Classification = parse_factor(Classification, levels = cOrderAlt)) %>%
mutate(gene = parse_factor(gene, levels = rev(gene.order))) %>%
mutate(type = parse_factor(
type,
levels = c(
"Wnt ligand", "Wnt target", "Other stemness factors",
"Stem-adjacent\nsecretory cell factors"
))) %>%
ggplot() +
geom_point(
aes(Classification, gene, size = pct, colour = scaled.mean.exp)
) +
facet_grid(type ~ tissue, space = "free", scales = "free") +
scale_color_gradient(low = "white", high = "darkred") +
scale.func(range = c(0, 12), limits = c(scale.min, scale.max)) +
guides(
size = guide_legend(
title = "% expressed", title.position = "top", title.hjust = 0.5
),
colour = guide_colorbar(
title = "Scaled mean expression", title.position = "top",
title.hjust = 0.5, barwidth = 10
)
) +
theme_bw() +
theme(
axis.text.x = element_text(angle = 90, hjust = 1),
axis.title = element_blank(),
legend.position = "top",
legend.justification = "center"
)## Joining, by = "Sample"p
ggsave(
plot = p,
filename = '../figures/MGA.enge.stemness.dotplot.pdf',
device = cairo_pdf,
height = 250,
width = 250,
units = "mm"
)Relative frequency of cell types in singlets vs. deconvoluted multiplets
singlets <- c(table(getData(cObjSng, "classification")))
singlets <- singlets / sum(singlets)
deconv <- colSums(adjustFractions(cObjSng, cObjMul, sObj))
deconv <- deconv[match(names(singlets), names(deconv))]
deconv <- deconv / sum(deconv)
if(!identical(names(singlets), names(deconv))) stop("name mismatch")
p <- tibble(
class = names(singlets),
singlet.freq = singlets,
multiplet.freq = deconv
) %>%
ggplot() +
geom_point(aes(singlet.freq, multiplet.freq, colour = class), size = 3) +
scale_colour_manual(values = palette('total')[order(names(palette('total')))]) +
xlim(min(c(deconv, singlets)), max(c(deconv, singlets))) +
ylim(min(c(deconv, singlets)), max(c(deconv, singlets))) +
geom_abline(slope = 1, intercept = 0, lty = 3, colour = "grey") +
labs(x = "Singlet relative frequency", y = "Multiplet relative frequency") +
guides(colour = guide_legend(title = "Cell Type")) +
theme_few()
p
ggsave(
plot = p,
filename = '../figures/MGA.enge.sngMulRelFreq.pdf',
device = cairo_pdf,
height = 160,
width = 250,
units = "mm"
)plotSwarmCircos(
sObj, cObjSng, cObjMul, weightCut = 0, classOrder = classOrder.MGA('total'),
classColour = palette('total')[classOrder.MGA('total')], h.ratio = 0.85
)## Joining, by = "class"
Detected duplicates - quadruplicates.
Only ERCC estimated cell number max 4.
Weight cutoff = 10.
# adj <- adjustFractions(cObjSng, cObjMul, sObj, binary = TRUE)
# samples <- rownames(adj)
# rs <- rowSums(adj)
# keep <- rs == 2 | rs == 3 | rs == 4
plotSwarmCircos(
sObj, cObjSng, cObjMul, weightCut = 10,
classOrder = classOrder.MGA('total'), classColour = palette('total')[classOrder.MGA('total')], h.ratio = 0.85,
alpha = 1e-3
)## Joining, by = "class"
pdf('../figures/MGA.enge.circos.pdf', width = 12, height = 10, onefile = FALSE)
plotSwarmCircos(
sObj, cObjSng, cObjMul, weightCut = 10,
classOrder = classOrder.MGA('total'), classColour = palette('total')[classOrder.MGA('total')], h.ratio = 0.85,
alpha = 1e-3
)## Joining, by = "class"invisible(dev.off())sessionInfo()## R version 3.6.1 (2019-07-05)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] patchwork_0.0.1 ggthemes_4.2.0 forcats_0.4.0
## [4] stringr_1.4.0 dplyr_0.8.3 purrr_0.3.2
## [7] readr_1.3.1 tidyr_0.8.3 tibble_2.1.3
## [10] ggplot2_3.2.1 tidyverse_1.2.1 CIMseq.data_0.0.1.4
## [13] CIMseq.testing_0.0.2 CIMseq_0.3.0.2
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-141 matrixStats_0.55.0 lubridate_1.7.4
## [4] RColorBrewer_1.1-2 httr_1.4.1 gmodels_2.18.1
## [7] tools_3.6.1 backports_1.1.4 R6_2.4.0
## [10] lazyeval_0.2.2 BiocGenerics_0.30.0 colorspace_1.4-1
## [13] withr_2.1.2 tidyselect_0.2.5 gridExtra_2.3
## [16] compiler_3.6.1 cli_1.1.0 rvest_0.3.4
## [19] xml2_1.2.2 labeling_0.3 scales_1.0.0
## [22] digest_0.6.20 rmarkdown_1.15 pkgconfig_2.0.2
## [25] htmltools_0.3.6 rlang_0.4.0 GlobalOptions_0.1.0
## [28] readxl_1.3.1 rstudioapi_0.10 shape_1.4.4
## [31] farver_1.1.0 generics_0.0.2 jsonlite_1.6
## [34] gtools_3.8.1 magrittr_1.5 Rcpp_1.0.2
## [37] munsell_0.5.0 S4Vectors_0.22.1 viridis_0.5.1
## [40] stringi_1.4.3 EngeMetadata_0.1.2 yaml_2.2.0
## [43] ggraph_2.0.0 MASS_7.3-51.4 plyr_1.8.4
## [46] Rtsne_0.15 grid_3.6.1 parallel_3.6.1
## [49] gdata_2.18.0 listenv_0.7.0 ggrepel_0.8.1
## [52] crayon_1.3.4 lattice_0.20-38 haven_2.1.1
## [55] graphlayouts_0.5.0 circlize_0.4.8 hms_0.5.1
## [58] zeallot_0.1.0 knitr_1.24 pillar_1.4.2
## [61] igraph_1.2.4.1 pso_1.0.3 reshape2_1.4.3
## [64] future.apply_1.3.0 codetools_0.2-16 stats4_3.6.1
## [67] glue_1.3.1 evaluate_0.14 modelr_0.1.5
## [70] vctrs_0.2.0 tweenr_1.0.1 cellranger_1.1.0
## [73] gtable_0.3.0 RANN_2.6.1 polyclip_1.10-0
## [76] future_1.14.0 assertthat_0.2.1 xfun_0.9
## [79] gridBase_0.4-7 ggforce_0.3.1 broom_0.5.2
## [82] tidygraph_1.1.2 googledrive_1.0.0 viridisLite_0.3.0
## [85] globals_0.12.4